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I.    INTRODUCTION 

This  thesis  models  numerically,  the  effects  of  a  pressure  disturbance  incident  on 
a  fluid  loaded  reinforced  plate  using  finite  difference  methods.  The  program  produced 
explores  the  reaction  of  a  plane  wave  striking  the  plate  (hull  of  a  ship)  at  different 
frequencies,  for  different  plate  lengths  and  for  various  spacings  of  the  reinforcing  beams 
(or  members). 

Figure  1  shows  the  problem  being  modeled.  The  figure  shows  a  cross  section  of 
a  reinforced  flat  plate  with  water  on  the  top  side  of  it  and  a  vacuum  below  it.  For  this 
problem,  it  is  assumed  that  the  plate  is  infinite  in  the  x  direction  along  the  horizontal 
axis.  Since  the  plate  structure  is  periodic  in  space,  the  model  will  only  look  at  one 
section  of  this  infinite  plate.  The  behavior  at  other  sections  is  assumed  to  be  the  same. 
We  will  use  three  fundamental  equations  in  the  model.  The  first  equation  is  the  wave 
equation,  the  second  is  the  linear-inviscid-force  equation,  and  the  third  is  the  equation 
of  motion  for  a  thin  plate.  The  coupling  of  these  equations  at  the  fluid-plate  interface 
makes  the  problem  very  difficult  to  solve  analytically. 

Chapter  II  derives  the  relationships  among  these  three  equations  as  they  are  used 
in  the  model.  Assumptions  used  in  the  program  are  spelled  out  in  this  chapter.  Chapter 
II  also  includes  the  steps  necessary  to  scale  each  equation  to  facilitate  its  use  in  the 
computer  program. 


Chapter  III  explains  the  operation  of  the  FORTRAN  computer  program  used  to 
model  the  problem.  The  program  takes  as  an  input  the  incident  and  specularly-reflected 
sound  waves  and  determines  the  scattered  pressure  due  to  the  fluid-structure  coupling. 
To  make  best  use  of  the  Naval  Postgraduate  School's  mainframe  computer  the  main 
program  is  divided  into  a  number  of  smaller  subroutines.  These  are  also  described  in 
Chapter  III. 

Chapter  IV  discusses  the  results  of  the  program  runs.  It  also  reviews  the 
weaknesses  in  the  model  and  model  limitations. 
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Figure  1.    Cross  Section  of  Hull 


II.    THEORY 

A.      THE  WAVE  EQUATION 

It  is  known  that  sound  travels  through  seawater,  or  any  fluid  medium,  as  a  pressure 
wave.  This  is  possible  because  seawater  is  a  compressible  fluid  which  is  affected  by 
pressure  disturbances.  To  approximate  the  movement  of  sound  in  a  fluid  medium  we 
have  used  the  linearized,  lossless  wave  equation  for  pressure  as  derived  in  Kinsler,  Frey, 
Coppens  and  Sanders  [Ref.  l:p.  104], 


il„T. 


1     *PT  ™ 


c)    dt2 


where  cf  is  the  fluid  sound  speed  and  pT  is  the  total  pressure. 

The  pressure  wave  can  be  split  up  into  various  parts,  each  of  which  is  a  solution 
of  the  linear  wave  equation.  In  this  case  we  want  to  be  able  to  solve  for  pressures  that 
would  be  present  in  the  fluid  which  modify  the  pressure  resulting  from  specular  reflection 
of  an  incident  plane  wave  from  an  acoustically  hard  surface.  This  known  specularly 
reflected  pressure  wave  is  combined  with  the  computed  scattered  pressure  wave  to  make 
up  the  total  reflected  pressure  wave. 

Symbolically,  we  have: 


pT=p,+pR+ps,  (2-2) 


where  pT  is  the  total  pressure,  p1  is  the  incident  pressure,  pR  is  the  specularly  reflected 
pressure,  and  ps  is  the  scattered  pressure.  The  total  radiated  pressure  is  a  superposition 
of  the  reflected  and  scattered  pressures. 

To  eliminate  the  specularly  reflected  wave  from  the  equations  of  motion  we  apply 
the  boundary  condition 

dp^  +  dp^l     =Q  (2.3) 

In  this  case  the  normal  derivatives  of  the  incident  and  specularly  reflected 
waves  must  cancel  at  the  surface  of  the  plate.  The  vertical  coordinate  Y  is  equal  to  zero 
at  the  surface. 

We  assume  that  the  source  of  the  incident  pressure  is  far  from  the  surface  so  that 
a  plane  wave  approximation  is  adequate  to  model  the  incident  pressure,  p1. 

The  time  harmonic  incident  and  specularly  reflected  plane  waves  are  of  the  form: 


!»(£=— 0 


pI(x,y,f)  =  e       C/  (2.4) 


and 


icj( -t) 


pR(x,y,t)  =  e       Cf 


(2.5) 


where 


and 


dl- 


COS0; 

-sin0, 


<?*= 


COS0; 

sin0, 


(2.6) 


(2.7) 


are  the  unit  direction  vectors  of  the  two  plane  waves. 

Both  the  incident  wave,  Equation  2.4  and  the  reflected  wave,  Equation  2.5  must 
satisfy  Equation  2.1.  Therefore,  because  of  superposition,  ps  must  also  satisfy  the  wave 
equation, 


c]  at2 


(2.8) 


B.      THE  LINEAR  INVISCID  FORCE  EQUATION 

The  linear  inviscid  force  equation,  taken  from  Kinsler,  Frey,  Coppens  and  Sanders 
[Ref.  l:p.  104],  is  basically  the  inviscid  Navier  Stokes  equation  in  which  nonlinear  terms 


have  been  neglected.  It  is  valid  for  acoustic  disturbances  of  small  amplitude,  and  can  be 
written 

pM*-*'.  (2'9) 

}  dt 

where  pf  is  the  constant  equilibrium  density  of  the  fluid,  and  u  is  the  velocity  of  the 
fluid. 

Since  we  assume  no  cavitation  of  the  fluid  or  penetration  of  the  fluid  at  the  fluid- 
plate  interface,  the  component  of  the  fluid  velocity  normal  to  the  interface  and  the  plate 
transverse  velocity  are  equal  at  the  fluid-plate  interface  at  all  times  t.  Letting  w  represent 
the  transverse  displacement  of  the  plate  (y  direction),  we  have 

p>=_&L     at  y=Q  (210) 

}  dt2         dy 
Through  use  of  Equations  2.2  and  2.3  the  total  pressure  pT  can  be  replaced  by  the 
scattered  pressure  ps: 

ppL.JeL         at    v=0.  (2.11) 

fdt2         dy 

C.      EQUATION  OF  MOTION  FOR  A  THIN  PLATE 

To  properly  understand  what  is  happening  on  and  within  the  thin  plate,  a  reference 
system  must  be  assigned.  In  Figure  2  we  see  that  the  vertical  plane  will  be  designated 
the  y-axis  in  which  the  thin  plate  experiences  transverse  displacement  w.  The  horizontal 
plane  will  be  the  x-axis  which  has  longitudinal  displacement  u.  To  set  up  our  equations, 


we  assume  that  the  centerline  (axis)  of  the  beam  stays  perpendicular  to  its  cross  section 
during  bending. 
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Figure  2  Plate  Reference  System 

As  put  forth  by  Kinsler,  Frey,  Coppens  and  Sanders  [Ref.  l:p.  58],  Hooke's  law 
states  that  if  the  strain  is  small  the  stress  is  proportional  to  it.  The  stress  on  the  plate  is 
described  by  the  equation 


°«  =  -y 


e   a2 


w 


1-v2  dx2 


(2.12) 


where: 


^  is  the  stress  on  the  plate,  E  is  the  Young's  modulus,  and  v  is  the  Poisson's  ratio. 


Consider  a  small  section  of  the  plate  of  area  AxAz  with  moments  and  forces  applied 
as  shown  in  Figure  3: 
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Figure  3.    Forces  and  Moments 
Note  that  there  are  no  forces  or  moments  applied  which  have  z  dependence.  The 
motion  therefore  is  assumed  to  be  two-dimensional.    Since  the  plate  is  not  rotating,  a 
force  balance  is  performed  on  the  moments  and  moment  arm,  neglecting  terms  of  order 


Ax2,  giving  us 

™-=F,  (2.13) 

dx 

where: 

M  =  bending  moment  on  the  element, 

q  =  load  on  the  plate  element  (uniform  for  a  normally  incident  wave), 

and  F  =  flexural  force  applied  to  the  element. 

Equating  forces  in  the  y-direction  applied  to  the  section  of  plate  yields 

ps—  (hAxAz)  =  —  AxAz  -  qAxAz,  (2. 14) 

dt2  d* 

where: 

S  =  cross  sectional  area  of  the  plate, 

pg  =  density  of  beam  section, 

h  =  plate  thickness. 

Dividing  out  the  AxAz  terms  yields 


-q+—  =  psh — _  (2.15) 


dx    '5    dt2 
Using  the  definition  of  the  moment 


\  1-v2  dx2 


2 

where  the  moment  of  inertia  is  found  from 


(2.16) 
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h 


I=fy2dy  =  !Lt  (2.17) 

j,  12 


we  obtain  the  equation  for  deflection  of  a  plate 


EI    d*w       ,  d*w  ,0  1Qv 

+  ph — -  =  -q.  (2.18) 


l-v2ac4      s    dt2 

The  loading  of  the  thin  plate  is  due  to  pressures  exerted  by  a  fluid  in  contact  with 
the  thin  plate.    Rewriting  Equation  2.18  including  the  pressure  terms  yields: 

8V      u#w      _T  (2-!9) 


D^?^ 


in  which 


D=      Eh*      ,  (2.20) 

12(l-v2) 

is  the  flexural  rigidity  and  h  is  the  plate  thickness,  Junger  and  Feit  [Ref.  2:p.  194]. 


D.      THE  SCALED  EQUATIONS 

1.      The  Scaling  Factors 

In  order  to  use  the  three  governing  Equations  2.1,  2.9  and  2.19  in  the  finite 
difference  code,  they  must  be  scaled  or  nondimensionalized.  The  steps  are  as  follows. 
The  independent  variables  x  and  t  are  scaled  through  use  of  the  fundamental  length  L 
and  the  wave  frequency  w,  respectively. 


10 


x  =  -t      and       y  =  2,  (2.21) 

L  L 

where  x  and  y    are  the  scaled  spatial  coordinates,  and  L  is  the  one  half  the  distance 
between  stiffeners  measured  in  the  x-direction. 
The  equation  for  the  scaled  time  is 


(2.22) 


(2.23) 


where  u>  is  the  frequency  of  excitation  of  the  incident  plane  wave. 
We  then  have 

_8__    a    _i  a 

dx    d(Lx)    L  dx 
and 


dt=      J    _G)  ft'  (2.24) 

d(— ) 

Additionally,  the  dependent  variables  w  and  p  must  be  scaled  as  follows.  For 
the  displacement  of  the  plate  (in  the  y-direction)  we  put 

w=^,  (2.25) 

W 

where  w    is  the  scaled  vertical  displacement  of  the  plate  and  W    is  the  displacement 
scaling  factor.    For  pressure,  put 
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Lu2pfW     P0 


where  p  is  the  scaled  pressure  and  P0  =  LurpfW  is  the  pressure  scaling  factor. 

The  wave  numbers  for  waves  propagating  in  the  fluid,  and  along  the  plate  are 
given  as  follows, 


*,=-,  (2.27) 


v 
V 


and 


co2p  h  \  (2.28) 

K=[— rH4, 


r 


D 


respectively. 

Additional  parameters  resulting  from  the  scaling  analysis  are 

q.£.jl.  <2-29> 

k;    ">c 

where  ft  is  the  squared  ratio  of  fluid  to  plate  wave  numbers  and  foG,  the  coincidence 
frequency  of  the  plate,  is  defined  by 

G  D 

We  also  define  the  dimensionless  quantity 
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e=-^L,  (2.31) 

where  e  is  a  dimensionless  quality  which  is  a  measure  of  the  fluid  loading  on  the  plate. 

2.      The  Wave  Equation 

The  wave  equation  governing  the  scattered  fluid  pressure  was  found  to  be: 


^Ps=\—y.  (2.32) 


9 


2dt2 


Scaling  the  above  equation  yields: 


°,-Sn-jiF*ntT2P*f  (2-33) 


which  reduces  to  the  desired  non-dimensional  form, 


Pz+P^^fpTT^PTT-  (2-34) 


cf 


3.      The  Linear  Inviscid  Force  Equation 

The  linear  inviscid  force  equation 

cPw  _  _  dps 
9/~d^~  "dy 

in  scaled  form  reduces  to 


(2.36) 


^77-Py-  (2-37) 
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4.      Thin  Plate  Equation 

The  dimensional  equation  of  flexural  rigidity  is: 


d*w       ,  cPw 


By  substituting  in  the  scaling  factors  it  becomes: 

Dy^----+ph  co2L4wivj--=  -p  TL5u2pjw,  (2.39) 

or,  after  more  simplification,  it  becomes 


-+k<L*w-----L5ekti(B+2B\  (2-40) 

Q 


This  is  the  scaled  form  of  the  thin  plate  equation  that  will  be  used  in  the 
computer  program.   Note  that  at  y  =  0,  pR  equals  p1.   Hence,  pT  has  been  replaced  with 
P  +  2p'. 

5.      Stiffeners 

In  the  present  analysis  the  stiffeners  will  be  treated  in  an  idealized  way.   The 
end  conditions  we  will  assume  are  imposed  at  the  point  where  the  stiffeners  are  located 
will  be  equivalent  to  those  of  a  clamped  plate  for  which  the  boundary  conditions  are 
w  =  0  at   x  =  ±  1 ,  and  wx  =  0   x  =  ±  1 . 

E.      THE  FLUID  DOMAIN 

As  stated  before,  it  is  assumed  the  plate  on  which  the  sound  wave  impinges  is 
infinite.  That  assumption  makes  it  possible  to  establish  the  following  conditions.  The 
boundary  conditions  on  the  section  of  plate  we  are  looking  at  will  be  periodic.    That  is 

14 


to  say,  for  normal  incidence  of  p\  displacement  on  the  plate  at  x  =  -1  will  be  the  same 
as  at  x  =  1.  For  arbitrary  angles  of  incidence,  the  pressures  at  x  =  ±1  differ  only  by 
a  scaling  factor  determined  by  the  angle  of  incidence.  Figure  4  and  the  following 
equations  will  show  how  the  time  harmonic  solution  is  dealt  with  in  the  problem. 
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Figure  4.    Side  View  of  Fluid  Domain 

To  understand  how  the  fluid  pressure  is  disturbed,  the  wave  equation  is  considered 
for  which  a  prescribed  pressure  is  applied  from  the  fluid  to  the  surface  of  the  plate.  This 
uncoupled  problem  can  be  solved  analytically  by  the  method  of  separation  of  variables. 

The  scaled  wave  equation  is 
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P    +P    =  a2Ptt.  (2.41) 

xx        yy  tt 

Applying  periodic  boundary  conditions  results  in 

p(-l$)=p(l$)      and     px(-l$)=px(l$).  (2.42) 

The  forced  pressure  at  the  plate  surface  is  given  by  the  boundary  conditions 

/>(*A0  =£(*)*-""         where      g(-\)=g(l)      and      gx(-l)=gx(l),       <2-43) 
where: 
g(x)  =  prescribed  pressure  applied  to  the  fluid  by  the  plate. 

Assuming  a  time  harmonic  solution  with  equivalent  frequency  dependence  to  that 
of  the  forcing  function,  the  scattered  pressure  can  be  expressed  as 

p(JJ,0=4)(iF,y)e-/r.  <2-44) 

By  substitution  of  Equation  2.44  into  Equation  2.41   we  obtain  the  reduced 
Hehmholtz  equation  for  4>(\,y), 

4>-+(j)--+a24>  =  0,  (2-45) 

with  boundary  conditions: 

<K-l,)T)  =  4>(l,y)      and      (|>x(-l,y)=<Mlj)  (2-46) 

and 
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<K*,0)  =  s(x).  (2-47> 

Henceforth  we  will  be  dropping  the  bars   from  all  equations,   giving  us  the 
following, 

^-^./aVO,  (2-48) 

<K-ioO=<Kij),  (2-49) 

<M-ioO=4>x(ioO.  (2-5°) 

♦,(-i.y)-Wiy).  (2-51) 

To  solve  for  0  the  standard  method  of  "Separation  of  Variables"  is  applied  to  the 
equation.   First  we  write  0  as  a  product  of  functions  of  one  variable  each, 

$=X(x)Y(y).  (2-52) 

Upon  substitution  into  the  Helmholtz  equation  and  separation  of  variables,  we  find  a 
solution 


X{x)=C  einnx.  (2.53) 

When  solved  for  Y,  two  answers  result: 

YJy)=Aeiy^-^      when      a>nn,  <2-54) 


and 
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Yn(y)=Ae-y^n2"2-al      when        a<rm  .  (2-55) 

These  solutions  are  chosen  to  avoid  the  physically  implausible  solution  that  the  scattered 
pressure  as  y  approaches  infinity  will  blow  up  represent  an  incoming  wave. 
Combining  our  solutions  we  have 

N 

*=5Z*«(X'^=5Z  aneb,nx*i^a2-n2"2)y+  £    a^fanx-)VnV-«2  (2.56) 

Note  that  as  y  approaches  infinity  the  decaying  terms  become  exponentially  small.   Now 
from  our  boundary  condition  at  y  =  0,    o>n  can  be  found  such  that 

i  ' 
g(x)=^aneinnx        with        an=-  fg(x)e-innxdx.  (2-57) 

2  J 


n=0 


1 


It  is  important  to  note  that  in  the  numerical  solution  of  the  fully  coupled  problem, 
the  radiation  waves  must  be  properly  modeled  as  y  approaches  infinity.  The  numerical 
domain  cannot  be  semi-infinite  in  extent  in  the  y  direction.  It  must  be  truncated  at  some 
value  of  y,  which  we  call  yB,  where  a  boundary  condition  must  be  imposed  to  ensure 
proper  modeling  of  the  radiated  pressure. 

F.       THE  PERIODICITY  PROBLEM 

To  ensure  that  a  wave  with  any  angle  of  incidence  can  be  correctly  modeled  by  the 
computer  program  the  periodicity  must  be  worked  out.  If  we  were  to  have  a  wave 
normally  incident  on  the  plate  we  would  expect  a  solution  for  the  pressure  which  is 
periodic  in  x  to  be 
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*=E*n(y)tf""  (2-58) 

To  account  for  an  arbitrary  angle  of  incidence,  Equation  2.58  needs  to  be  modified 
in  the  following  manner.  A  correction  factor  which  takes  into  account  the  phase  change 
at  x  =  -1  vice  x  =  1  must  be  added.  The  modified  equation  for  the  scattered  pressure 
is, 


*=E*nO0e'"nI}e~'*C"6/X,  <2-59> 


where  0,  is  the  angle  of  propagation  of  the  incident  pressure  wave. 
Continuing  with  this  argument  one  can  write 

Y^-kcosdj+rnz.  (2.60) 

Substitution  of  this  x  dependence  into  the  Helmholtz  equation  yields 

^+(*2-yX=<>.  (2-61) 

dyl 

Substituting  these  relationships  back  into  our  equation  for  <f>  gives  us 

N  -A/-1 

♦  ■  E  v*"*"^  E  +E  >«/'"'''.  (2-62> 


n  =  -M  -oo         N*\ 


where 
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P.- 


fi^n      *>Y„  (263) 

fin1*     *<Y„ 


and  in  which  the  first  sum  includes  all  values  of  n  for  which  k  >  7n  and  the  latter  two 
sums  account  for  the  remaining  modes. 

G.      FINITE-DIFFERENCE  MODELING 

1.      The  Taylor  Series 

Refer  to  Numerical  Solution  of  Partial  Differential  Equations:  Finite 
Difference  Methods.  Chapter  I,  by  G.  D.  Smith  [Ref.  3]  for  a  more  complete  explanation 
of  finite  difference  methods. 

All  three  governing  equations  will  be  center  differenced  with  truncation  errors 
of  order  (Ax2,  At2). 

Finite  difference  approximations  for  ux,  uxx,  uXJUCX  are 


w<— ^w2).  (264) 

2/1 


(uJi~  ,     "l+(Xh2)  (2.65) 


and 


(u.^^-%,^)^^  p.^ 


where  u4  =  u(Xi). 
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2.      Differencing  of  the  Wave  Equation 

We  have  already  shown  that  the  wave  equation  is  of  the  following  form 

P     +P     =    k?L2Ptt,  (2.67) 

xx  yy  J  tt' 

and  that  we  have  the  following  conditions, 

x,  =  iAx,  yj  =  jAy,  tn  =  nAt,  (2.68) 

where  i,  j,  and  n  are  integers.    This  allows  us  to  describe  the  pressure  as  follows, 

P"j  =  P(xityptn)         and         Ax  =  Ay=h.  (2-69) 

So  by  substituting  into  Equation  2.66  we  arrive  at  the  final  form, 


3.      Finite  Differencing  of  the  Linear  Inviscid  Force  Equation 

We  have  already  shown  that  the  linear  inviscid  force  equation  is  of  the 
following  form 

dy 

This  equation  is  applied  at  y  =  0  or  when  j  =  0  in  the  numerical  domain. 
The  displacement  w  does  not  depend  upon  y,  so  we  introduce  W"  =  W(x„  tn )  and  write 
the  differenced  form  of  the  equation  as: 


21 


W^  =  -At2  [  *J±J±±\+2W?-W*-\  (2-72> 

2Ay 

This  equation  is  solved  for  the  Pn, A  which  is  exterior  to  the  numerical  domain. 
It  can  then  be  substituted  into  the  differenced  form  of  the  wave  equation  to  solve  for 
Pn+1i0  ,  the  pressure  on  the  plate  fluid  interface  at  time  level  n  +  1  solving  for  ViA  we 
have 


p£_x  =  1±1[W"*X-2W?  +  W-~X]+  P.V  (2.73) 


Ar2 

Note  that  Wn+1i  is  first  needed  to  solve  for  Pn; .,  which  means  a  prescribed  order  must  be 
followed  at  the  fluid  plate  interface.  First  Wn+1i  is  found,  the  P"^.,,  and  finally  Pn+1ii0. 

4.      Finite  Differencing  the  Plate  Equation 

The  last  of  the  major  equations  that  we  must  put  in  the  form  of  finite 
differences  is  the  thin  plate  equation.    The  scaled  thin  plate  equation  is: 


W^ktVW-  -  l£^(p+2P')  .  <2-74) 


XXXX        P  I  I  f-v 

Differencing  this  equation  results  in  an  expression  for  the  displacement  W  at 
time  level  n+ 1  in  terms  of  the  preceding  two  time  levels  and  the  incident  plane  scattered 
pressure  at  time  level  n. 
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w>  ~-—t —      — 

~p 


k-L  "  (2.75) 


-^^-ip^ip  %o,tn)]  ^w'-wr1 


5.      Finite  Differencing  Side  Boundaries 

When  differencing  the  wave  equation  at  the  side  boundaries  for  the  pressure 
we  have  a  similar  problem  to  that  found  in  applying  the  wave  equation  at  the  fluid  plate 
interface.  The  problem  is  that  there  are  values  of  the  pressure  needed  which  are  exterior 
to  the  numerical  domain.  Periodicity  will  be  used  as  an  additional  condition  to  eliminate 
this  problem.  We  have  the  scattered  pressure  in  the  form: 


f,,M£<t>„(w)ei'""k~'hcose'  <2-76> 


What  we  are  doing  here  is  taking  the  pressure  at  the  node  i  =  M  +  l,  x  =  Ax 
and  replacing  the  pressure  i  =  -M  +  l  and  x  =  -1  +Ax.   From  Equation  2.76  we  obtain 

K.^E^yX"'"""^^'^'  (2-77) 

K.u-^'^T  *&*-&*-*>-*•  (2-78) 

and  we  also  have: 

'WE  <t>„(v>*""""4xV  'm'"'")cae■  <2-79> 

n 

Now  we  can  make  the  following  two  statements: 
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K>\re  '*Axoos%  -,kcosQ'J2  ♦,(-i)"«*,"A?  (2-80> 

Combining  the  previous  statements  we  get: 

n  -2«fcose,  p„  (2.82) 

An  equivalent  argument  used  at  the  left  boundary  yields: 

n  likcosQ      n  (2.83) 

r-M-\J    e  rM-\J  V  ' 

6.      Finite  Differencing  the  Radiated  Boundary. 

The  final  edge  of  the  numerical  domain  is  the  artificial  boundary  which  must 
allow  free  passage  of  radiated  pressure  disturbances.  It  must  be  far  enough  away  so  that 
the  evanescent  modes  can  be  neglected.  To  see  how  the  boundary  condition  is  derived 
we  start  with  the  scattered  pressure  ps  and  try  to  derive  an  analytical  expression  for  the 
pressures  normal  direction  at  the  boundary,  as  is  done  in  Scandrett  and  Kriegsmann 
[Ref.  4].  This  type  of  radiation  boundary  condition  is  often  referred  to  as  a  non-local 
radiation  boundarv  condition. 

J 

From  previously  we  have 

N 

ps(x,y,t)  =e  ~u  {  ]T  ane'Pn>'+Y"x+  evanesent  modes),  (2.84) 

n—M 

where 
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v   =-ikcosQ,+nn, 

'  n  I  ' 


(2.85) 


P."R 


(2.86) 


and 


VpS=kWr 


(2.87) 


Additionally  in  our  problem  we  have 


k=k,L     and       k/= — , 
cf 


(2.88) 


but  as  y  approaches  infinity  ps  is  equivalent  to  only  the  radiated  modes.   We  include  the 
time  dependence  into  the  a„  coefficient  such  that, 


A' 


PS~Y,  an(t)e 


•(P-VY.r*) 


n  =  -M 


and  in  the  following  steps  use  orthogonality  to  find  an(t)ei^ny.    We  have 


(2.89) 


fe -,Y"> s{x,y,t)dx  =  £  «„(0e'P"7« '^e^dx  (2-90) 


n=M 


-1 


and 


ix[-i£cos0,+rt7t  -(-iicosOi+mii)] 
_  u[nit-mn] 


(2.91) 


We  can  therefore  write 


where  |  =  dummy  variable  of  integration. 
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"„(')«  "■'  "|  /«"""'/' s(«W)d«.  (2-92> 


2-, 


Noting  the  similarity  to  zero  in  the  following  equation 


y  [JL  +  Q  ±}a  «  V,"y*v)  =  0.  (2-93) 


Which  can  be  rewritten  as 


^r4  E  P,/'WHI^({W«-  0.  (2.94) 


Remember,  these  equations  pertain  to  the  artificial  boundary  which  is  modeling  infinity. 
For  this  problem  y  is  equal  to  yB  at  the  artificial  boundary. 

When  we  start  to  finite  difference  Equation  2.94,  the  artificial  boundary  and 
use  the  trapezoidal  rule  of  integration  we  obtain: 

pIR  pM  ft  IX 

,,.         *TT  E  P.*  E  V"-'"""^1-^'))-  =0  (2-95) 

in  which  we  take: 

1 


l  =  ±IX 

*r\  2 

1      l*±IX. 


(2.96) 


Rewriting  Equation  2.95  we  have 
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«"  Dm  ,       IX  N 


^y1L^Li+JL^6/[^  ^/^^W^-P^yO  (2-97) 


2/j  4Ar/=n     B=_A, 


or 


p  *     _  p w         /x 

,,     '£^S'-tf)-'  (2-98) 


Where  the  square  matrix  A  is  defined  as 

A  =l!Lby  pV^'l  (2-99) 

Combining  the  above  equations  gives  us: 

i 

pm        _pm  y,      .     (pm*l_pm-U_Q  (2.100) 


/x 

£ 

From  the  differenced  form  of  the  wave  equation  applied  at  j   =   IY  (the  artificial 
boundary)  we  have: 

p%  -  -K^-^n^K^p^K^ ^~^   <2- I01> 

h  k  h  kz 

Combining  Equations  2.100  and  2.101  and  eliminating  Pni,iY+i  we  obtain: 
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IX  2  (2-102^ 

ehhhc  lJ  IJY    •         k2h2  ,JY 

Now  we  need  to  put  this  equation  into  a  form  which  makes  it  possible  for  us 
to  see  the  matrixes  needed  to  solve  the  problem.    First  we  write: 


ri,IY 

+ 
h2k 

,   £   AiflJY 
l=-IX 

= 

My 

Af2 
h2k2 

IX 
1=-IX 

f ijy 

hV-l-,jr  ■ 

1 

-IJY' 

¥lPn 

(2-4 

Ar2      . 

define 

i  J1  fr  ,=/ 


and 


Equation  2.102  can  then  be  written  in  matrix  form  as 


where 


(2.103) 


6.,  =  <I    I   lr\  (2-104) 


Au-»u+— An  (2.105) 


B.^-b.^-^-A..  (2.106) 


AP;+1=BP^  +  C  (2-107) 


28 


C<  =  ■^[p.".,n-+P<-W+2^-i]+(2-4-^i)P(,"„,  (2.108) 

/tat  h  kl 

Then  the  equation  for  the  pressure  vector  at  the  next  time  level  is  given  by: 

P£1=A-1[BPn-l+C}  (2-109) 

The  pressure  at  Pn+ViY  is  tnen  assigned  by  taking  the  i*  component  of  the 
vector  P0"1"1^.  In  solving  the  matrix  Equation  2.107  we  use  the  LINPACK  library 
Gaussian  elimination  programs  (CGESO  and  CGESL). 
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III.    DESCRIPTION  OF  PROGRAM 

A.  GENERAL 

The  computer  program  written  for  this  thesis  simulates  the  effects  of  a  pressure 
disturbance  on  a  fluid  loaded  reinforced  thin  plate.  A  detailed  description  of  the  theory 
behind  the  program  can  be  found  in  Chapter  II.  A  pressure  wave  in  a  fluid  medium  is 
moving  towards  a  reinforced  plate  which  is  flat,  and  is  supported  by  stiffeners  which  are 
at  a  constant  spacing  on  the  non-fluid  side  of  the  plate.  The  non-fluid  side  is  assumed 
to  be  a  vacuum.  When  the  pressure  sound  wave  hits  the  plate  there  will  be  some  (a  very 
small  amount)  distortion  of  the  plate  in  the  vertical  direction.  This  displacement  will 
cause  a  disturbance  in  the  fluid  medium.  The  scalded  pressure  amplitude  for  different 
frequencies  and  beam  spacing  will  then  be  plotted  and  the  magnitude  of  the  propagating 
modes  will  be  recorded. 

The  program  stores  the  time  dependent  wave  equation  with  time  harmonic  forcing 
by  allowing  transient  effects  to  radiate  away  from  the  numerical  domain.  When  a  steady; 
state  has  ben  achieved  the  code  automatically  stops  the  calculations. 

B.  SUBROUTINES 

To  make  effective  use  of  the  Postgraduate  School's  mainframe  computer,  the 
program  has  been  broken  down  into  a  number  of  smaller  programs  (i.e.  subroutines) 
which  can  be  more  efficiently  programmed  and  run. 
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The  main  program  is  used  to  coordinate  the  subroutines  and  to  pass  information 
from  one  subroutine  to  another.    The  following  are  descriptions  of  the  subroutines. 

1.  Subroutine  INIT 

The  INIT  subroutine  initializes  the  computer  program.  In  this  subroutine  we 
set  up  our  two  main  matrices.  The  first  of  these  is  the  "P"  matrix,  which  represents  the 
fluid  pressure  in  the  fluid  domain.  The  second  is  the  "W"  matrix,  this  represents  the 
displacement  of  the  plate  surface  which  bounds  the  fluid  medium. 

Although  we  have  the  ability  to  change  any  number  of  parameters  in  the 
program,  we  will  concentrate  on  two.  These  are  the  beam  spacing,  and  the  frequency 
of  the  pressure  wave  impinging  on  the  plate.  Both  of  these  values  are  assigned  in  this 
subroutine.  The  INIT  subroutine  sets  all  of  our  counters  and  subscripts  to  zero  at  this 
time  and  initializes  the  matrices  which  will  be  needed  to  model  the  artificial  boundary. 

2.  Subroutine  DOMAIN 

The  DOMAIN  subroutine  calculates  the  pressure  at  all  nodes  in  the  interior 
of  the  fluid  medium  at  a  subsequent  time  level  in  terms  of  the  preceding  pressure  values. 
The  subroutine  calculates  the  pressure  at  the  next  time  step  by  use  of  the  differenced 
wave  Equation  2.70. 

The  equation  above  must  be  modified  to  model  what  is  happening  at  the  side 
boundaries  of  the  fluid  domain.  As  was  stated  earlier,  our  model  is  looking  at  a  small 
section  of  an  "infinite"  fluid  plate  interface.  This  allows  us  to  use  the  pressure  in  the 
next  section  to  model  what  is  happening  in  the  section  in  which  we  are  interested.   This 
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is  necessary  since  otherwise  we  would  have  to  use  "undefined"  points  in  our  differenced 
wave  equation  at  the  edges  of  the  fluid  domain. 

3.  Subroutine  FLPL 

The  FLPL  subroutine  models  the  fluid-plate  interface  along  the  top  of  the 
reinforced  thin  plate.  Again,  with  finite  differencing  and  the  assumption  that  we  are  only 
looking  at  a  small  section  of  an  infinite  interface,  we  are  able  to  model  the  effects  of  the 
pressure  wave  along  the  entire  fluid-plate  interface.  This  subroutine  is  also  where  we 
introduce  our  driving  function  into  the  computer  program.  This  function  is  what  starts 
things  moving  towards  a  steady  state.  Although  this  is  not  the  most  complicated 
subroutine  in  the  entire  program  it  is  the  one  which  has  proved  to  be  the  hardest  to  make 
work  as  predicted. 

4.  Subroutine  ARTF 

The  ARTF  subroutine  is  the  most  complicated  of  any  subroutine  in  the 
computer  program.  It  models  the  artificial  boundary  of  the  fluid  medium  at  infinity. 
Again,  this  is  accomplished  with  finite  differencing  and  matrix  operations. 

5.  Subroutine  ALPHA 

The  ALPHA  subroutine  calculates  the  magnitude  of  the  radiating  modes  of  the 
pressure  wave.    This  data  is  gathered  in  the  results  section  of  this  study. 

6.  Subroutines  SHIF  and  FINN 

The  SHIF  subroutine  moves  the  entire  problem  to  the  next  time  level.  This 
is  done  after  the  other  subroutines  have  completed  all  the  required  calculations  needed 
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to  satisfy  our  domain  and  all  of  our  boundary  conditions.  After  exiting  the  SHIF 
subroutine  the  entire  program  is  run  again  with  the  exception  of  the  INIT  subroutine. 
The  FINN  subroutine  is  a  check  to  make  sure  that  we  have  reached  a  steady  state 
situation  before  we  take  any  measurements  of  the  plate  displacement  or  of  any  fluid 
disturbances.  Figure  5  is  a  graphical  representation  of  the  displacement  convergence 
factor  approaching  zero.  This  was  arrived  at  with  a  assigned  frequency  of  3500  Hz  and 
L  =  1.0. 

7.      Subroutine  CLOSE 

The  CLOSE  subroutine  is  used  to  send  the  data  generated  in  other  parts  of  the 
main  program  to  the  plotting  devices. 

C.      PROGRAM  VERIFICATION 

Verification  of  the  program  was  accomplished  by  the  following  method.  We  solved 
the  waveguide  problem  (section  2E)  for  a  known  value  of  g(x).  We  let  g(x)  =  sin  7rx  and 
eixx.  This  was  accomplished  by  replacing  the  FLPT  subroutine  with  a  subroutine  named 
TEST.  As  before  the  program  was  run  until  a  steady  state  solution  was  reached.  These 
results  were  checked  against  analytical  solutions  and  found  to  compare  very  well.  For 
g(x)  =  sin  xx  we  graphically  reproduced  the  expected  results.  For  g(x)  =  elxx  we 
calculated  the  coefficients  an  and  found  a,  —  1  and  all  other  an's  ~  0.  The  small  error 
found  in  the  calculation  can  be  attributed  to  truncation  errors  in  the  finite  difference 
approximations.   These  results  match  with  analytical  calculations. 
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CONVERGENCE  GRAPH 
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Figure  5.    Convergence  Graph 
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IV  RESULTS  AND  CONCLUSIONS 
A.      INITIAL  CONDITIONS 

To  properly  apply  the  computer  program  generated,  inputted  data,  or  conditions, 
must  be  within  the  tolerances  of  the  main  equations  used.  This  is  demonstrated  by  the 
use  of  .01  meters  as  the  thickness  of  the  thin  plate  and  a  frequency  range  from  500  Hz 
to  4000  Hz.  Other  inputted  data,  including  the  densities  of  steel  and  water  were  taken 
directly  from  Kinsler,  Frey,  Coppens  and  Sanders  [Ref.  l:pp.  461-463]. 

The  choice  of  the  Ax  and  Ay  terms  is  based  on  the  need  to  have  a  significant 
number  of  data  points  in  the  waveguide  domain  that  will  give  an  accurate  picture  of  what 
is  happening.  This  is  also  true  for  the  time  step  chosen,  we  need  to  make  sure  that  the 
At  was  small  enough  to  keep  the  problem  from  blowing  up.  A  Courant-Lewy-Friedrichs 
stability  condition  must  be  satisfied  for  the  domain  which  limits  the  size  of  At. 
Additionally  as  was  mentioned  before,  the  computer  program  was  run  for  the  same 
amount  of  time  for  each  trial,  this  time  period  is  long  enough  to  ensure  that  a  steady  state 
solution  was  being  achieved. 
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B.      PROGRAM  OUTPUT  DATA 

The  data  obtained  from  the  computer  program  demonstrates  the  following  major 
points.  When  the  frequency  is  increased  you  get  more  propagating  modes  from  each 
particular  plate  length.  When  L  is  changed  by  a  specific  amount  you  do  see  a  difference 
in  the  magnitude  and  shape  of  the  deformation  graphs. 

Our  first  table  shows  the  calculated  magnitudes  of  the  propagating  modes  of  the 
pressure  wave  at  an  assigned  frequency  of  750  Hz  at  different  beam  spacings. 
Table  I.    MAGNITUDE  OF  RADIATING  MODES. 


For  a  frequency  of  750  Hz 

with  L  =  0.8      a0  =  .467808127 

with  L  =  1.0      a0  =  .484098911 
a,  =  .307484090 

with  L  =  1.2      a0  =  .232444227 
a,  =  .367214561 


In  addition  to  the  data  in  Table  I,  we  have  in  Figures  6  through  8  a  graphical 
representation  of  the  scalded  pressure  amplitude.  The  pressure  was  measured  at  the 
artificial  boundary  ,  yB  at  a  frequency  of  750  Hz.  The  horizontal  axis  is  of  a  length  of 
2L  for  each  of  the  graphs  starting  at  x  =  -1  and  going  to  x  =  1.  The  numbering  on  the 
horizontal  axis  represents  each  of  the  41  pressure  nodes  used  in  the  finite  difference 
code. 
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Figure  6.   Sample  Pressure  Graph 
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L  =  1.0        FREQUENCY  =  750  HZ 
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Figure  7.   Sample  Pressure  Graph 
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L  =  1.2        FREQUENCY  =  750  HZ 
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Figure  8.    Sample  Pressure  Graph 
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Table  II  is  a  list  of  the  magnitudes  of  the  an's  for  program  runs  with  a  frequency 
of  3500  Hz. 

Table  II.    MAGNITUDE  OF  RADIATING  MODES. 


For  a  frequency  of  3500  Hz 

withL  =  0.8  a0  =  .451181054 
a,  =  .323620617 
a2  =  .088650584 

withL  =  1.0  a0  =  .365310729 
a,  =  .260265648 
a2  =  .062102031 
a3  =  .060141488 
a,  =  .076183497 

withL  =  1.2       c*o  =  .252027571 

a,  =  .165709317 

a,  =  .047427550 

a3  =  .046469088 

aA  =  .049672558 

as  =  .043821129 

a6  =  .067380726 


From  the  data  above  we  observe  that  at  a  frequency  of  3500  Hz  there  is  a 
substantial  change  in  the  magnitude  of  the  propagating  modes  as  L  is  increased  from  0.8 
to  1.2.  It  can  be  seen  from  other  runs  that  the  higher  the  frequency  the  larger  the 
differences  in  the  magnitudes  of  the  propagating  modes. 

Figures  8  through  10  show  the  scaled  pressure  amplitude  data  obtained  from  the 
program  runs  with  an  assigned  frequency  of  3500  Hz.  It  is  observed  that  changing  of 
the  length  L  on  which  the  Ax  is  scaled  has  an  effect  on  not  only  the  amplitude  of  the 
pressure  wave  but  also  their  shape. 
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Figure  9.   Sample  Pressure  Graph 
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L  =  1.0        FREQUENCY  =  3500  HZ 
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Figure  10.   Sample  Pressure  Graph 
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L  =  1.2       FREQUENCY  =  3500  HZ 
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Figure  11.   Sample  Pressure  Graph 
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Table  III  is  a  list  smaple  energy  calculations  taken  at  three  frequencies  inputted  into 
the  computer  program.  The  data  presented  demonstrates  that  for  theses  frequencies  and 
lengths  run  in  the  computer  program  that  the  fundamental  mode  contains  the  largest 
amount  of  energy.  Also  evident  is  that  the  energy  contained  in  the  fundamental  mode 
starts  becomes  more  constant  as  the  frequency  is  increased  for  different  lengths  of  L. 
This  data  is  in  agreement  with  the  theory  of  sound  propagation. 
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Table  III   ENERGY  CALCULATIONS 


Percent  of  Energy  in 
Propagating  Mode 

For  a  frequency  of  750  Hz 

with  L  =  0.8        E0  =  100% 

withL  =  1.0        E0  =    99% 
E,  =     1% 


Percent  of  Energy  in 
Propagating  Mode 

For  a  frequency  of  3500  Hz 


with  L  =  0.8 


E0  =  66% 
E,  =  32% 
E9  =     2% 


with  L  =  1.2 


E0  =    77% 
E,  =    23% 


with  L  =  1.0 


For  a  frequency  of  1750  Hz 


with  L  =  0.8 


with  L  =  1.0 


with  L  =  1.2 


E0  =  77% 

E,  =  23% 

E0  =  68% 

E,  =  31% 

Ej  =  1% 

E0  =  64% 

E,  =  32% 

Ej  =  3% 

E3  =  1% 


with  L  =  1.2 


E0  = 

65% 

E,  = 

32% 

£2  = 

2% 

E3  = 

1% 

E4  = 

* 

E0  = 

64% 

E,  = 

27% 

E2  = 

2% 

E3  = 

2% 

E4  = 

2% 

E5  = 

1% 

E6  = 

2% 

*  denotes  much  less  then  1  % 
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Table  IV  shows  the  magnitude  of  the  primary  propagating  mode  for  each  of  the 
different  frequencies  and  lengths  inputted  into  the  computer  program. 
Table  TV   MAGNITUDES  OF  PRIMARY  MODES 

Magnitude  of  a0's 


L  =  0.8 


L  =  1.0 


L  =  1.2 


Frequency 
750  Hz 

1000  Hz 

1250  Hz 

1500  Hz 

1750  Hz 

2000  Hz 

2250  Hz 

2500  Hz 

2750  Hz 

3000  Hz 

3250  Hz 

3500  Hz 

3750  Hz 


.467808127   .484098911   .367214561 


.6182562721   .444052100   .262355387 


.553467214 


.486700475 


.42002557; 


.457091272   .343211770 


.498957574   .415586054   .320898831 


.404151559   .431289613   .338643909 


.433078885   .316784620 


.394520938   .400619507   .252146780 


.385663450   .392564356   .267684340 


.385663450   .392564354   .267684340 


,353338659   .234915197 


.369634628   .226448834   .194508076 


.451181054   .365310729   .252027571 


.446322470   .330691218   .226533830 
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The  data  obtained  from  the  computer  program  and  represented  in  Table  III  shows 
that  there  is  always  a  difference  in  the  magnitude  of  the  fundamental  propagating  mode 
for  changes  in  L  at  all  frequencies  comprised  within  the  computer  runs.  The  difference 
in  the  magnitude  of  the  primary  propagating  mode  between  the  runs  with  a  length  "L" 
of  1.2  and  the  runs  with  a  length  L  of  1  are  greater  then  the  difference  between  the  runs 
with  lengths  of  0.8  and  1.0.  In  the  frequency  range  covered  when  L  is  increased  from 
1.0  to  1.2  there  is  always  a  drop  in  the  magnitude  of  the  fundamental  propagating  mode. 
This  is  not  always  true  when  L  is  decreased  to  0.8  from  1.0.  Although  the  magnitude 
of  these  increases  are  smaller  then  the  magnitude  of  the  decreases  this  demonstrates  that 
by  increasing  the  spacing  of  the  reinforcing  beams  will  not  always  result  in  a  decrease 
in  the  amplitude  of  the  fundamental  propagating  mode  in  this  frequency  range.  From  this 
result  we  can  conclude  that  additional  factors  play  an  important  role  in  determining  the 
magnitude  of  the  propagating  modes. 

Because  the  model  is  allowed  to  run  until  a  steady  state  solution  is  reached  there 
is  no  appreciable  difference  in  the  absolute  pressure  measured  at  the  artificial  boundary 
for  any  of  the  frequencies  observed.  However,  it  is  evident  that  as  the  length  of  the  plate 
"L"  is  increased  there  is  a  slight  drop  in  the  pressure  measurement.  This  is  true  for  all 
frequencies.  Figure  12  is  a  sample  plate  deformation  graph  for  L  =  1.0  and  a  frequency 
of  1000  Hz. 

In  conclusion,  the  computer  program  produces  results  as  expected  for  the  different 
frequencies  of  the  plane  wave  and  lengths  of  the  thin  plate.  The  superpositioning  of  the 
reflected  plane  wave  has  showed  that  as  L  is  increased  there  is  a  larger  number  of 
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propagating  modes  accompanied  by  a  change  in  the  amplitude  of  each  mode.  Additional 
results  could  be  explored  for  different  angles  of  deflection  as  well  as  different  plate 
thicknesses.  This  will  require  modification  of  the  computer  program  as  it  stands  to 
gather  data  for  additional  angles  and  thicknesses. 
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Figure  12.    Deformation  Graph 
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APPENDIX:  COMPUTER  PROGRAM 


FILE:  280CT     FORTRAN   Al 
C 

C     .  THIS  PROGRAM  SIMULATES  THE  DISTRUBANCE  IN  A  FLUID 

C       AFTER  A  SOUND  WAVE  HAS  HIT  A  PLATE. 

C 

C      P  =  PRESSURE  (EXCESS)       D  =  FLEXURAL  RIGIDITY 

C      X  =  HORIZONTIAL  POSITION    Y  =  VERTICAL  POSITION 

C      W  =  DEFORMATION  KF  =  FLUID  WAVE  NUMBER 

C    SPD  =  WAVE  SPEED  KP  =  PLATE  WAVE  NUMBER 

C   FREQ  =  FREQUENCY  OMEGA  =  RATIO  OF  (KF/KP)**2 

C 

C      COUNTERS i   I   J   N 

C 

C         THE  VALUES  IN  THE  PARAMETER  STATEMENTS  ARE  CONSTANTS 

PARAMETER  (IX  =  20,  IY  =  80) 

PARAMETER  (H  =  .025,  DELTAT  =  .0325   ) 

PARAMETER  (NDIM  =  2XIX+1) 
C  THE  FOLLOWING  ARE  COMPLEX  ARRAYS 

COMPLEXX16    P(  -IX:  IX,  Oi  IY,  3),  W(  -IX.  IX,  3) 

COMPLEX       LB,  FACR 

COMPLEX       A(  -IX:  IX,  -IX:  IX) 

COMPLEX       AA(  NDIM,  NDIM) 

COMPLEX       BB(  NDIM,  NDIM) 

COMPLEX      C(  -IX:  IX) 

COMPLEX       Z(  NDIM  ) 

C0MPLEXX16    Fl 

COMMON/SS/P,  W 

COMMON/TT/KP,KF,INFA,  THETAI,L,  X,  PI, 
1  EISP,  OMEGA,  T 

COMMON/CC/LB,  FACR 

COMMON/EE/WGTWO,  WGTHE,  WGTOT,  CON 

COMMON/PT/PGTWO,  PGTHE,  PGTOT 

COMMON/HH/MINN,  MAXX 

COMMON/II/AA,  BB 

COMMON/KK/IPVT 

COMMON/MM/NN 

COMMON/OO/C 

COMMON/FF/A 

COMMON/UP/AL 

COMMON/RN/RUNS 

COMMON/FEL/U,V,TP 

COMMON/SN/FREQ,  SPD 

COMMON/GF/WGGRF, PGGRF 
C 

INTEGER   RUNS,  IPVT(NDIM),  SAM 

PARAMETER  (SAM  =  1000) 

REAL      KF,  KP,  INFA,  EISP,  OMEGA,  L,PI 

REALX8    WGTWO, WGTHE, WGTOT, CON 

REALX8    PGTWO, PGTHE, PGTOT 

REAL      WGGRF(  SAM  ),  PGGRF(  SAM  ) 

REAL      THETAI,  FREQ,  SPD 

REAL      X,  T 

PI  =  3.141592654 

CON   =  .015D0 
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c 
c 

c 
c 


WGTHO  =  0.0D0 

WGTHE  =  O.ODO 

PGTWO  =  O.ODO 

PGTHE  =  O.ODO 

c 

c 

CALL  XUFLOW  (0) 
CALL  INIT 
RUNS  =  0 

20 

CALL  FLPL 
CALL  DOMAIN 

RUNS  =  RUNS  +  1 

CALL  ARTF 

IF  (RUNS  .GT.  1999)  THEN 

GOTO  60 

ELSE 

GOTO  61 

ENDIF 

60 

CALL  FINN 

61 

T  =  T  +  DELTAT 
CALL  SHIF 

IF  (RUNS  .EQ.  3000)  THEN 
GOTO  80 
ELSE 
GOTO  20 
ENDIF 

80 

CALL  ALPHA 
CALL  CLOSE 

WRITE(2,767)THETAI,KF 

767  FORMATdX, 'THETAI    =    ',F12.6,3X,'      KF    =    SF12.6) 
WRITE(2,768)0MEGA,EISP 

768  FORMATdX, »    OMEGA    =    ' ,  F12  .6,  3X,  »  EISP    =    ',F12.7) 
WRITE(2,196)DELTAT,    KP 

196    FORMATdX, 'DELTAT    =    ',F12.6,3X,'      KP    =    ',F12.6) 

WRITE(2,487)H,L 
487    FORMATdX,'  H    =    »,F12.6,3X,*         L    =    ',F6.3) 

WRITE(2,735)FREQ,SPD 
735    FORMATdX, 'FREQ    =    ' ,  F12  .6  ,  3X,  'SOUND   SPEED   =    SF12.6) 

WRITE(2,629)PI 
629    FORMATdX,'      PI    =    »,F12.7) 

WRITE(2,919) 
919    FORMATdX,  » FACR    =    •) 

WRITE(2,*)    FACR 
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c 
c 
c 
c 
c 

WRITE(2,222) 

222  FORMATUX,  'FINAL  DATA' ,//,  IX,  ■        X     ', 
1   '     W      »,  '      P     ') 

DO  333  I  =  -IX,  IX,  1 

XX    =    IXH 
333        WRITE(2,666)XX,CDABS(W(I,2)),CDABS(P(I,IY,3)) 
666         FORMATUX, 5F10. 5) 

HRITE(2,223)RUNS 

223  <      FORMATUX, 'RUNS    =    ',17) 

WRITE(2,678) 
678         FORMATUX, 'BELOW    IS   THE    FINAL    WGTOT,    AND   PGTOT') 
WRITE(2,X)WGT0T 
WRITE(2,x)PGT0T 

STOP 

END 
C 

SUBROUTINE  INIT 
C 

C         THIS  SETS  UP  THE  INITIAL  CONDITIONS  OF  THE 
C         PROBLEM  INCLUDING  THE  DOMAIN  AND  THE 
C         ARTIFICIAL  BOUNDARY 
C 
C         THE  VALUES  IN  THE  PARAMETER  STATEMENTS  ARE  CONSTANTS 

PARAMETER  (IX  =  20,  IY  =  80) 

PARAMETER  (H  =  .025,  DELTAT  =  .0325) 

PARAMETER  (NDIM  =  2XIX+1) 
C  THE  FOLLOWING  ARE  COMPLEX  ARRAYS 

COMPLEXX16    P(  -IX:  IX,  0:  IY,  3),  W(  -IX*  IX,  3) 

COMPLEX       LB,  FACR 

COMPLEX       A(  -IX:IX,  -IX:  IX) 

COMPLEX       AA(  NDIM    ,  NDIM    ),  BB(  NDIM    ,  NDIM 

COMPLEX       Z(   NDIM   ) 

COMMON/SS/P,  W 

COMMON/TT/KP,KF,INFA,  THETAI,L,  X,  PI, 


FILE:  280CT     FORTRAN   Al 


1  EISP,  OMEGA,  T 

COMMON/CC/LB,  FACR 
COMMON/FF/A 

COMMON/GG/GAMMA,  BETTA 
COMMON/HH/MINN,  MAXX 
COMMON/II/AA,  BB 
COMMON/KK/IPVT 
COMMON/FEL/U,V,TP 
COMMON/SN/FREQ,  SPD 

REAL  THETAI,  INFA,  KF,  EISP,  OMEGA,  TIN,  L,  KP 

REAL  ART,  K,  XI,  XIL,  RCOND 

REAL  GAMMA(-40:  40),  BETTA  (-40:40),  PI 

REAL  SPD,  FREQ 

REAL  X,  T 

INTEGER  I,  J,  N,  LJ,  NN,  U,  V,  TP 

INTEGER  NMODMIN,  NMODMAX,  MINN,  MAXX 

INTEGER  IPVT(NDIM) 
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u 

TP  =  41 

LB  =  CMPLX(0.0,  1.0) 

FACR  =  CMPLXd.O,  0.0) 

SPD  =  1500. ODO 

FREQ  =  3750. ODO 

N  =  0 

1  =  0 

J  =  0 

NN  =  1 

L  =   1.2D0 

THETAI  =   PI/2. ODO 

T  =  DELTAT 

NMODMIN  =  0 

NMODMAX  =  0 

KF  =  (2D0XPIXFREQ*L)/SPD 

c 

\ 

GT  =  ((195000000000. Ox( . 05**3) )/12(l-( .28XX2)) 

c 

KP  =  ((((FREQ*2*PI)X*2)X7700XH)/FR)XX.25 
KP  =  40.7534975147D0 
OMEGA  =  (KF/KP)x«2 
INFA  =  (DELTATX*2)/((KFX*2)X(LXX2)X(HXX2)) 
DO  12  I  =  -IX,IX,1 
DO  11  J  =  0,IY,1 

P(I,J,1)  =  CMPLX(0.0,  0.0) 
P(I,J,2)  =  CMPLX(0.0,  0.0) 
PCI, J, 3)  =  CMPLX(0.0,  0.0) 

11 

CONTINUE 

c 

12 

CONTINUE 

DO  21  I  =  -IX,  IX,  1 

W(I,1)  =  CMPLX(0.0,  0.0) 

H(I,2)  =  CMPLX(0.0,  0.0) 

c 

W(I,3)  =  CMPLX(0.0,  0.0) 

c 
c 
c 

21 

CONTINUE 

THIS  SETS  UP  THE  ARTIFICAL  BOUNDRAY 

15 


K  =  KF*L 

TIN  =  -KXCOS(THETAI) 

DO  15  N  =  -40,  40,  1 

GAMMA(N)  =  TIN  +  PIXN 
IF  (  ABS(GAMMA(N)) 
NMODMIN  =  MIN(N, 
NMODMAX  =  MAX(N, 
ENDIF 
CONTINUE 

MINN  =  NMODMIN 
MAXX  =  NMODMAX 


LE.  K)  THEN 

NMODMIN) 

NMODMAX) 


DO  18  N  =  MINN,  MAXX,  1 

BETTA(N)  =   SQRTUKXX2)  -  (GAMMA(  N)x*2) ) 
18  CONTINUE 

DO  78  I  =  -IX,  IX,  1 
DO  79  LJ  =   -IX,  IX,  1 

A(I,LJ)  =  CMPLX(0. 0,0.0) 
79    CONTINUE 
78  CONTINUE 

ART  =  (2X(HXX2))/(4XDELTAT) 
DO  300   I  =  -IX,  IX,  1 

DO  301   LJ  =  -IX,  IX,  1 
XI  =  IXH 
XI L  =  L  JXH 

IF  (ABS(LJ)  .EQ.  IX)  THEN 
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DEL  =  .5 
ELSE 
DEL  =  1 
ENDIF 
DO  302  N  =  MINN,  MAXX,  1 
A(I,LJ)  =  A(I,LJ)  +  ART*DEL 
1  XBETTA(N)*CEXP(LBXGAMMA(N)X(XI-XIL)) 

302  CONTINUE 
301  CONTINUE 
300  CONTINUE 


DO 
DO 


IF  (I 


420 
421 


U  =  1, 
V  =  1, 
I  = 
LJ  = 
AA(U,V)  = 
BB(U,V) 
EQ.  LJ) 
AA(U,V) 
BB(U,V) 


NDIM,  1 

NDIM,  1 

-IX+U-1 
=  -IX+V-1 
=   INFAX(A(I,LJ)) 
=   INFAX(A(I,LJ)) 
THEN 

=  AA(U,V)  +  1.0 
=  BB(U,V)  -  1.0 


ENDIF 
421  CONTINUE 
420  CONTINUE 

CALL  CGECO 

RETURN 

END 


(AA, NDIM, NDIM,  IPVT,  RCOND,Z) 


SUBROUTINE  FLPL 

THIS  SUBROUTINE  DEFINES  THE  FLUID-PLATE  INTERFACE 
THIS  IS  WHERE  Y=0  ALONG  THE  TOP  OF  THE  PLATE. 

THE  VALUES  IN  THE  PARAMETER  STATEMENTS  ARE  CONSTANTS 
PARAMETER  (IX  =  20,  IY  =  80) 
PARAMETER  (H  =  .025,  DELTAT  =  .0325) 
PARAMETER  (NDIM  =  2*IX+1) 

THE  FOLLOWING  ARE  COMPLEX  ARRAYS 


IX:  IX,  0:  IY, 
FACR,  FTEMP 


C0MPLEXX16    P( 
COMPLEX       LB, 
C0MPLEXX16    Fl 
COMMON/SS/P,  W 

COMMON/TT/KP,KF,INFA,  THETAI,L, 
1  EISP,  OMEGA,  T 

COMuGN/'CC'' L  B ,  FACR 
COMMON/RN/RUNS 


3),  W(  -IX:  IX,  3) 


X,  PI, 


REAL     KP,  L,  OMEGA,  EISP,  KF,  INFA,  PI,  THETAI 

REAL     X,  T 

REALX4   SQ,SW,BW1,BW2,BW3,BW4 

INTEGER  I,  J,  RUNS 


EISP  =  0.13101864566D0 

J  =  0 

SQ  =  (DELTAT/(HXKF))XX2 

SW  =  (DELTAT/(HXKP)XX2)XX2 

BW1  =  2-6XSW 

BW2  =  4XSW 

BW3  =  EISPXDELTATXX2XKF/0MEGA 

BW4  =  -SW 

A  =  FL0AT(IX)/FL0AT(2XIX) 
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c  ~~ 

C         THIS  LOOP  COMPUTES  THE  INTERIOR  NODES 
C 

DO  10  I  =  -IX  +  2,  IX  -2,  1 

y  —     T  £  U 

WCI,3)  =  BW1XWCI,2)  -WCI,1)  +  BW2XCWCI+1,2)+WCI-1,2)) 

1  +BW4xCWCI+2,2)+WCI-2,2)) 

2  +BW3X(P(I,0,2)+  F1(X,T)) 

PTEMP      =  -(2*H/DELTAT*X2)K(U(I,3)-2*W(I,2)+W(I,1))  +PCI,1,2) 
PCI,J,3)  =  INFAXCPCI+1,J,2)+PCI-1,J,2)-4XPCI,J,2)+PCI,J+1,2)+ 
1  PTEMP     )+2XPCI,J,2)  -PCI,J,1) 

10  CONTINUE 
C 

C         THESE  LOOPS  COMPUTE  THE  NODES  AT  -19  AND  19 
C 

I  =  -IX+1 

Y  =  f — Ty+1 }^H 

WCI,3)  =  BW1*WCI,2)-WCI,1)  +BW2X(W(I+1,2)+WCI-1,2)) 

1  +BN«xCWCI+2,2)+WCI,2)) 

2  +BW3*CPCI,0,2)+  F1(X,T)) 

PTEMP     =  -C2xH/DELTATxx2)xCWCI,3)-2*WCI,2)+WCI,l))  +PCI,1,2) 
PCI, J, 3)  =  INFAXCPCI+1,J,2)+PCI-1,J,2)-4XPCI,J,2)+PCI,J+1,2)+ 
1  PTEMP      )+2xPCI,J,2)  -PCI,J,1) 

C 

I  =  IX  -1 

y  —     ( T Y_ 1 ) kU 

NCI, 3)  =  BW1*WCI,2)  -  HCl,l)  +  BW2XCWCI+1,2)+WCI-1,2)) 

1  +BU4XCWCI,2)+WCI-2,2)) 

2  +BW3xCPCI,0,2)+  F1CX,T)) 

PTEMP      =  -(2XH/DELTATXX2)X(W(I,3)-2*W(I,2)+H(I,1))  +PCI,1,2) 
PCI, J, 3)  =  INFAX(P(I+1,J,2)+P(I-1,J,2)-4*P(I,J,2)+P(I,J+1,2)+ 
1  PTEMP      )+2XPCI,J,2)  -PCI,J,1) 

C 

C  THESE  ARE  THE  END  POINTS 

C 

FACR=CEXP(2XLBXLXKF*C0S(THETAI)) 

I  =  -IX 

J  =  0 

W(I,3)  =  CMPLX(0. 0,0.0) 

PTEMP   =  PCI, 1,2) 

PCI, J, 3)  =  INFAX(P(I  +  1,J,2)  +  FACRXP(-I-1,J,2)-<4XP(I,J,2) 
1  +P(I,J+1,2)+PTEMP      )+2xP(I,J,2)  -P(I,J,1) 

FACR  =  CEXP(-2XLBXLXKFXC0S(THETAD) 

I    =    IX 

J    =    0 

WCI,3)    =    CMPLX(0. 0,0.0) 

PTEMP      =    P(I,1,2) 

PCI, J, 3)  =  INFAXCFACRxPC-I  +  l,J,2)  +  PCI-l,J,2)-<ixpci,  J, 2) 
1  +PCI,J+1,2)+PTEMP      )+2XPCI,J,2)  -PCI, J, 1) 

RETURN 

END 
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FILE:  280CT     FORTRAN   Al 

SUBROUTINE  DOMAIN 
C 
C        THE  VALUES  IN  THE  PARAMETER  STATEMENTS  ARE  CONSTANTS 

PARAMETER  (IX  =  20,  IY  =  80) 

PARAMETER  (H  =  .025,  DELTAT  =  .0325) 

PARAMETER  (NDIM  =  2XIX+1) 
C         THE  FOLLOWING  ARE  COMPLEX  ARRAYS 

C0MPLEXX16    P(  -IX:  IX,  0:  IY,  3),  W(  -IX:  IX,  3) 

COMPLEX       LB,  FACR 

COMMON/SS/P,  W 

COMMON/TT/KP,  KF,  INFA,  THETAI,  L,  X,  PI, 
1  EISP,  OMEGA  ,T 

COMMON/CC/LB,  FACR 
C 

REAL     KP,  KF,  INFA,  THETAI,  L,  PI,  EISP,  OMEGA 

REAL     X,  T 

INTEGER  I, J 
C 

DO  22  I  =  1-IX,  IX-1,1 
DO  21  J  =  1,  IY-1,  1 
P(I,J,3)  =  INFAX(P(I+1,J,2)+P(I-1,J,2)- 

1  4XP(I,J,2)+P(I,J+1,2)+P(I,J-1,2))+ 

2  2*P(I,J,2)-P(I,J,1) 

21  CONTINUE 

22  CONTINUE 

LB  =  CMPLX  (0.0,  1,0) 

LA  =  LXKFXCOS(THETAI) 

FACR  =  CEXP((2)*(LBXLA)) 

I  =  -IX 

DO  41  J  =  1,  IY-1   ,  1 

P(I,J,3)  =  INFAX(P(I+1,J,2)+  FACRxp(-I-l,J,2) 

1  -4XP(I,J,2)+P(I,J+1,2)+P(I,J-1,2))+ 

2  2XP(I,J,2)-P(I,J,1) 
41  CONTINUE 

C 

LB  =  CMPLX  (0.0,  1.0) 

RA  =  LXKFXCOS(THETAI) 

FACR  =  CEXP((-2)X(LBXRA)) 

I  =  IX 

DO  31  J  =  1,  IY-1,1 

P(I,J,3)  =  INFAX(FACRXP(-I+1,J,2)+  P(I-1,J,2) 

1  -4XP(I,J,2)+P(I,J+1,2)+P(I,J-1,2))+ 

2  2XP(I,J,2)-P(I,J,1) 
31  CONTINUE 

RETURN 
END 
C 

r 
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SUBROUTINE  ARTF 

C         THE  VALUES  IN  THE  PARAMETER  STATEMENTS  ARE  CONSTANTS 

PARAMETER  (IX  =  20,  IY  =  80) 
•PARAMETER  (H  =  .025,  DELTAT  =  .0325) 

PARAMETER  (NDIM  =  2*IX+1)      ,MB%/- 
C  THE  FOLLOWING  ARE  COMPLEX  ARRAYS 

C0MPLEXX16    PC  "IXt  IX,  0:  IY,  3),  WC  -IX:  IX,  3) 

COMPLEX       LB,  FACR 

COMPLEX       AA(  NDIM   ',  NDIM    ),  BBC  NDIM    ,  NDIM 
COMPLEX       BCNDIM) 
COMPLEX       C(  -IX:  IX) 
COMPLEX       Z(   NDIM   ) 
COMMON/SS/P,  W 

COMMON/TT/KP,KF,INFA,  THETAI,L,X,  PI, 
1  EISP,  OMEGA  ,T 

COMMON/CC/LB,  FACR 
COMMON/ FF/ A 
COMMON/II/AA,  BB 

C 

SUBROUTINE  SHIF 
C 

C  THIS  MOVES  THE  PROGRAM  TO  THE  NEXT  TIME  LEVEL 

C 
C         THE  VALUES  IN  THE  PARAMETER  STATEMENTS  ARE  CONSTANTS 

PARAMETER  (IX  =  20,  IY  =  80) 

PARAMETER  (H  =  .025,  DELTAT  =  .0325) 

PARAMETER  (NDIM  =  2*IX+1) 
C  THE  FOLLOWING  ARE  COMPLEX  ARRAYS 

COMPLEXX16    P(  -IX:  IX,  0:  IY,  3),  W(  -IX:  IX,  3) 

COMPLEX       LB,  FACR 

COMMOH/SS/P,  W 

COMMON/TT/KP,KF,INFA,  THETAI,  L,  X,  PI, 
1  EISP,  OMEGA,  T 

COMMON/CC/LB,  FACR 
C 

REAL     KP,  L,  OMEGA,  EISP,  KF,  INFA,  PI,  THETAI 

REAL     X,  T 

INTEGER   I,  J 


III 
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COMMON/FEL/U,V,TP 

COMMON/KK/IPVT 

COMMON/OO/C 

REAL     KP,  L,  OMEGA,  EISP,  KF,  INFA,  PI,  THETAI 

REAL     X,  T 

INTEGER   IPVT(NDIM),  JOB  ,U,V,  TP 

TP  =  41 

JOB  =  0 

LB  =  CMPLX  (0.0,  1.0) 

LA  =  L*KF*COS(THETAI) 

RA  =  LXKFXCOS(THETAI) 

DO  42  I  =  -IX+1,  IX-1,  1 

C(I)  =  INFAX(P(I+1,IY,2)+P(I-1,IY,2)+2XP(I,IY-1,2)) 
1        +  (2  -  (4XINFA))XP(I,IY,2) 
42  CONTINUE 

I  =  -IX 

FACR  =  CEXP((2)X(LBXLA)) 

C(I)  =  INFAX(P(I+1,IY,2)+FACRXP(-I-1,IY,2)+2XP(I,IY-1,2)) 
1        +  (2  -  (4XINFA))XP(I,IY,2) 

I  =  IX 

FACR  =  CEXP((-2)X(LBXRA)) 

C(I)  =  INFAX(FACRXP(I-1,IY,2)+P(-I+1,IY,2)+2XP(I,IY-1,2)) 
1        +  (2  -  (4XINFA))XP(I,IY,2) 

DO  10  I  =    1,  NDIM,  1 

B(I)  =  INFAX(CMPLX(0.0,  0.0)) 
10  CONTINUE 

DO  20  I  =  1,NDIM,1 

DO  30  LJ  =  1,  NDIM,  1 

B(I)  =  B(I)  +  BB(I,LJ)XP(-IX+LJ-1,IY,1) 
30  CONTINUE 
20  CONTINUE 

DO  150  I  =    1,  NDIM,  1 

BCD  =  BCD  +  C(-IX+I-1) 
150  CONTINUE 

CALL  CGESL (AA, NDIM, NDIM, IPVT,B, JOB) 
DO  46  I  =  -IX,  IX,  1 

P(I,IY,3)  =  BdX+l  +  I) 
46  CONTINUE 
RETURN 
END 
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c 
c 


c 
c 


DO  19  I  =  -IX,  IX,  1 
DO  22  J  =  0,  IY,  1 
P(I,J,1)  =   P(I,J,2) 

22  CONTINUE 

19  CONTINUE 

DO  20  I  =  -IX,  IX,  1 
DO  23  J  =  0,  IY,  1 
PCI, J, 2)  =   P(I,J,3) 

23  CONTINUE 

20  CONTINUE 

DO  29  I  =  -IX,  IX,  1 
W(I,1)  =   W(I,2) 
29  CONTINUE 

DO  40  I  =  -IX,  IX,  1 
W(I,2)  =   W(I,3) 
40  CONTINUE 
RETURN 
END 

FUNCTION  FKXP,  TIME) 
PARAMETER  (IX  =  20,  IY  =  80) 
PARAMETER  (H  =  .025,  DELTAT  =  .0325) 
C0MPLEX*16    P(  -IX:  IX,  0:  IY,  3),  W(  -IX:  IX,  3) 
COMPLEXX16   Fl,  SI 
COMMON/SS/P,  W 
COMMON/RN/RUNS 

COMMON/TT/KP,KF,INFA,  THETAI,L,  X,  PI, 
1  EISP,  OMEGA,  T 

REAL     KP,  L,  OMEGA,  EISP,  KF,  INFA,  PI,  THETAI 
REAL     X,  T 
REAL*8  Q9,Q5,S 


IF(RUNS    .EQ.    1)    THEN 

S    =    0.0 

ELSE 
S      =    TIME+XP*KF*COS(THETAI) 

ENDIF 
Q9    =    0 
Q5    =    AMINKTIME,    18.0) 
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Q8  =  DEXP(2*Q5-6.0) 

Q7  =  Q8/(Q8+.1/Q8) 

SI  =  CMPLX(Q9,-S) 

Fl  =  Q7*CDEXP(S1) 

RETURN 

END 
C 
C 

.  SUBROUTINE  ALPHA 
C 

C        THIS  SETS  UP  THE  INITIAL  CONDITIONS  OF  THE 
C        PROBLEM  INCLUDING  THE  DOMAIN  AND  THE 
C        ARTIFICIAL  BOUNDARY 
C 
C        THE  VALUES  IN  THE  PARAMETER  STATEMENTS  ARE  CONSTANTS 

PARAMETER  (IX  =  20,  IY  =  80) 

PARAMETER  (H  =  .025,  DELTAT  =  .0325) 

PARAMETER  (NDIM  =  2*IX+1) 
C         THE  FOLLOWING  ARE  COMPLEX  ARRAYS 

COMPLEX*16    P(  -IX.  IX,  0:  IY,  3),  W(  -IX.  IX,  3) 

COMPLEX       LB,  FACR 

COMPLEX       A(  -IX:IX,  -IXi  IX) 

COMPLEX       AAC  NDIM    ,  NDIM    ),  BBC  NDIM    ,  NDIM 

COMPLEX       AL(  -40.  40) 


COMPLEX       LAND 
COMMON/SS/P,  W 

COMMON/TT/KP,KF,INFA,  THETAI,L,  X,  PI, 
1  EISP,  OMEGA,  T 

COMMON/CC/LB,  FACR 
COMMON/FF/A 

COMMON/GG/GAMMA,  BETTA 
COMMON/HH/MINN,  MAXX 
COMMON/II/AA,  BB 
COMMON/KK/IPVT 
COMMON/RN/RUNS 
COMMON/FEL/U,V,TP 
COMMON/SN/FREQ,  SPD 

REAL  THETAI,  INFA,  KF,  EISP,  OMEGA,  L,  KP 

REAL  XI,  X,  T 

REAL  GAMMA(-40.  40),  BETTA  (-40.40),  PI 

REAL  SPD,  FREQ 

INTEGER  I,  U,  V,  TP,  MINN,  MAXX,  IPVT(NDIM),  RUNS,  Y 

Y  =  IY/2 
WRITE(2,75) 
75   FORMATUX,'   ') 
HEN  =  RUNSXDELTAT 


C 


DO  10  M  =  MINN,  MAXX,1 

AL(M)  =  CMPLX(0.0,  0.0) 
DO  30  I  =  -IX,  IX,  1 
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XI  =  I*H 
IFC  ABS(I)  .EQ.  IX)  THEN 
DEL  =  .5 
ELSE 

DEL  =  1.0 
ENDIF 
LAND  =  CEXP(-LB*Y*HXBETTA(M)  +  (LB*HEN)) 

AL(M)  =  AL(M)  +  .5*LANDXHXDELXCEXP(-LB*GAMMA(M)XXI)*P(I,Y,3) 
30     CONTINUE 
10   CONTINUE 

WRITEC2,32) 
32   FORMATUX, 'ALPHA  nNnS  FOR  RADIATING  MODES') 
DO  111  M  =  MINN,  MAXX 
HRITE(2,X)CABS(AL(M)) 
111   CONTINUE 
RETURN 
END 
C 

SUBROUTINE  CLOSE 
C 

C         THE  VALUES  IN  THE  PARAMETER  STATEMENTS  ARE  CONSTANTS 
PARAMETER  (IX  =  20,  IY  =  80) 
PARAMETER  (H  =  .025,  DELTAT  =  .0325) 
PARAMETER  (XYX=  40) 
PARAMETER  (TNT=  1000) 
PARAMETER  (NDIM  =  2XIX+1) 
C  THE  FOLLOWING  ARE  COMPLEX  ARRAYS 

COMPLEXX16    P(  -IX:  IX,  0:  IY,  3),  W(  -IXt  IX,  3) 
COMPLEX       LB,  FACR 
C  THE  VALUES  IN  THE  COMMON  "BLOCKS"  ARE  VARIABLES 

COMMON/SS/P,  W 

COMMON/TT/KP,KF,INFA,  THETAI,L,  X,  PI, 
1  EISP,  OMEGA  ,T 

COMMON/CC/LB,  FACR 
COMMON/HH/MINN,  MAXX 
COMMON/GF/WGGRF, PGGRF 
C 

REAL      KP,  KF,  INFA,  THETAI,  L,  PI,  EISP,  OMEGA 

REAL      X,  T 

REAL      XX(  0.  41),  WPLOT(  -IX:  IX),  ZZ(0:1000) 
REAL      PPLOSC  -IX:  IX),  SMALL(IOOO) 
INTEGER   II,  RUNS,  SAM,  IR 
PARAMETER  (SAM  =  1000) 
REAL      WGGRF(  SAM  ),  PGGRF(  SAM  ) 
C 

DO  10  II  =  0,  XYX,  1 
XX(II)  =  II 

10  CONTINUE 
C 

DO  11  IR  =  0,  TNT,  1 
ZZ(IR)  =  IR 

11  CONTINUE 
C 

DO  30  I  =  1,1000,1 

SMALL(I)  =  WGGRF(I) 
30  CONTINUE 

DO  35  I  =  -IX,  IX,  1 

PPLOS(I)  =  CDABS(P(I,IY,3)) 
35  CONTINUE 

DO  45  I  =  -IX,  IX,  1 

WPLOT(I)  =  CDABS(W(I,3)) 
45  CONTINUE 
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c 
c 
c 
c 
c 

CALL  COMPRS 

CALL  PLOTD(ZZ, SMALL, 1000,  .TRUE.  , ' LINLIN' , »CONV  FACTOR*', 

1  'CONVERGENCE  GRAPHS', 'RUN  NUMBER  +  2000$', 

2  'CONVERGENCE  VALUES') 
C 

CALL  PL0TD(XX,PPL0S,41,  .TRUE.  ,' LINLIN' ,' FLUID  PRESSURES  * , 

1  'PI/2   LENGTH=1.2   FREQ=  3750$ ', 'HORIZONTAL  POSITIONS' 

2  'SCALED  PRESSURES') 

CALL  PLOTD(XX,WPLOT,41,  .TRUE.  ,' LINLIN ',» DISP  AT  3750  HZS ' , 

1  'SCALED  LENGTH  OF  1 . 2S ', 'SCALED  HORIZONTAL  POSITIONS', 

2  'SCALED  DISPLACEMENT  OF  PLATES') 
CALL  DONEPL 

RETURN 
,  END 
C 

SUBROUTINE  FINN 
C 

C        THIS  SUBROUTINE  MAKES  SURE  WE  HAVE  REACHED  THE  STEADY  STATE 
C         SOLUTION  OF  THE  PROBLEM,  IE  ' EOUALIBRIUM' 
C        THE  VALUES  IN  THE  PARAMETER  STATEMENTS  ARE  CONSTANTS 
C 

PARAMETER  (IX  =  20,  IY  =  80) 

PARAMETER  (H  =  .025,  DELTAT  =  .0325) 
C  THE  FOLLOWING  ARE  COMPLEX  ARRAYS 

C0MPLEXX16    p(  -IX:  IX,  0«  IY,  3),  W(  -IX«  IX,  3) 

COMMON/SS/P,  W 

COMMON/TT/KP,KF,INFA,  THETAI,  L,  X,  PI, 
1  EISP,  OMEGA  ,T 

COMMON/EE/WGTWO,  WGTHE,  WGTOT,  CON 

COMMON/PT/PGTWO,  PGTHE,  PGTOT 

COMMON/RN/RUNS 

COMMON/GF/WGGRF, PGGRF 
C 

REAL     KP,  L,  OMEGA,  EISP,  KF,  INFA,  PI,  THETAI 

REAL     X,  T 

REALX8   WGTWO, WGTHE, WGTOT, CON 

REALX8   PGTWO, PGTHE, PGTOT 

INTEGER  I,  RUNS,  SAM 

PARAMETER  (SAM  =  1000) 

REAL      WGGRF(  SAM  ),  PGGRF(  SAM  ) 
C 

DO   34  I  =   -IX+2,  IX-2,  1 

WGTWO   =    W(I,2)*DC0NJG(W(I,2))+  WGTWO 
WGTHE   =    W(I,3)*DCONJG(W(I,3))+  WGTHE 
WGTOT   =    (WGTHE  -  WGTWO)  /  WGTHE 
34  CONTINUE 
C 

DO   44  I  =   -IX   ,  IX   ,  1 

PGTWO   =    P(I,IY,2)XDC0NJG(P(I,IY,2))+  PGTWO 
PGTHE   =    P(I,IY,3)XDCONJG(P(I,IY,3))+  PGTHE 
PGTOT   =    (PGTHE  -  PGTWO)  /  PGTHE 
44  CONTINUE 
C 

IF  (RUNS.GT.  1999)  THEN 
WGGRF(RUNS-1999)  =  WGTOT 
PGGRF(RUNS-1999)  =  PGTOT 
ENDIF 
C 

C       WRITE(2,88) 

C   88         FORMATdX, 'WGGRF, PGGRF    IN    FINN') 
C  WRITE(2,X)    WGGRF(RUNS-1999) 

C  WRITE(2,X)    PGGRF(RUNS-1999) 

RETURN 
END 
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